Bayesian regression facilitates quantitative modelling of cell metabolism
This paper presents Maud, a command-line application that implements Bayesian statistical inference for kinetic models of biochemical metabolic reaction networks. Maud takes into account quantitative information from omics experiments and background knowledge, as well as structural information about kinetic mechanisms, regulatory interactions and enzyme knockouts. Our paper reviews the existing options in this area, presents a case study illustrating how Maud can be used to analyse a metabolic network and explains the biological, statistical and computational design decisions underpinning Maud.
Bayesian inference, kinetic models of cell metabolism, multiomics integration, ordinary differential equations, regulatory anlaysis
1 Introduction
A kinetic model of cellular metabolism aims to express what is known about a cellular process in the form of an in silico representation of the underlying network of chemical reactions. Kinetic models can be used to improve production of target molecules, determine regulatory networks (Christodoulou et al. 2018) and identify potential drug targets (DeBerardinis and Chandel 2016; Liberti et al. 2017). However, the use of kinetic models in practice is hindered by their dependence on noisy and uncertain information sources. Quantitative in vivo measurements of chemical abundances, and in vitro measurements relating to kinetic parameters, both contain vital information but are notoriously inaccurate (Siskos et al. 2017; Tabb et al. 2010; Lu et al. 2017). Practically useful kinetic modelling therefore requires a principled statistical approach that encompasses multiple possible model parameterisations.
Bayesian statistical inference can combine the structural information implicit in kinetic models with knowledge about metabolic parameters and information from omics measurements (Saa and Nielsen 2016; Gopalakrishnan, Dash, and Maranas 2020). However, kinetic models pose serious computational challenges for Bayesian inference (Gutenkunst et al. 2007; A. Raue et al. 2010).
The scope of a kinetic model is defined by a stoichiometric matrix, \(S\), in which rows represent metabolites, columns represent reactions, and matrix elements \(s_{ij}\) represent the stoichiometric coefficient of metabolite \(i\) in reaction \(j\). The change in metabolite concentrations is:
\[\begin{equation}\label{eq-1} \frac{dC}{dt} = S\cdot v - \mu\cdot C \end{equation}\]
where C represents a vector of metabolite concentrations, \(v\) is a vector of reaction rates, and \(\mu\) is the growth rate. The second term represents the dilution due to cell growth.
In a kinetic model, the rates, \(v\), are expressed as a function of the enzyme concentrations, \(E\), the metabolite concentrations, \(C\), and a set of parameters \(\theta\):
\[\begin{equation}\label{eq-2} v = f(C, E, \theta) \end{equation}\]
The parameters must include sufficient boundary concentrations and fluxes to solve \(\eqref{eq-1}\).
It is common to assume pseudo-steady state for metabolites, i.e., the rate of fluxes towards any metabolite is much greater than the rate of change in concentration, \(𝑣 \gg \frac{𝑑𝐶}{𝑑𝑡}\). Moreover, the dilution effect is assumed minimal, \(\mu\cdot C \ll \vec{v}\) (true unless the concentration is very high). Finally, the enzyme concentration is assumed constant for the period considered and hence part of the parameters.
Given these assumptions and a set of values for , a set of steady state metabolite concentrations and fluxes can be found by solving for \(C\) the algebraic equation:
\[\begin{equation}\label{eq-3} S\cdot f(C;\theta) = 0 \end{equation}\]
In a fermentation context, \(\eqref{eq-3}\) captures the rapid kinetics inside the cell, while another set of ODEs would be used to describe the external substrate and product concentrations, which could act as boundary parameters to \(\eqref{eq-3}\).
In the context of kinetic modelling, Bayesian inference is appealing because it allows uncertainty to be represented appropriately without sacrificing mechanistic accuracy. Measurement uncertainty can naturally be represented in a Bayesian measurement model, whereas the prior model can represent quantitative uncertainty about kinetic parameters. Finally, kinetic rate laws can be represented in Bayesian data generation models with arbitrarily high fidelity. See Gelman, Carlin, et al. (2020) for more details about Bayesian inference and Gelman, Vehtari, et al. (2020) for a discussion of practical Bayesian workflow.
Another advantage is that Bayesian inference problems are well-posed even when not all parameters are strongly identified. Sloppy parameter models in which measurable quantities are sensitive to combinations of parameters but not to individual marginal parameter values are ubiquitous in models of biological systems (Gutenkunst et al. 2007; White et al. 2016). The parameter correlation structure represents the set of potential models that describe the observed data. As we demonstrate later, capturing this correlation structure is difficult outside a fully Bayesian context.
Previous Bayesian kinetic models have either sacrificed mechanistic accuracy or have attempted to fit realistic kinetic models using obsolete or unreliable computational methods.
The most popular algorithm for fitting Bayesian statistical models is Markov Chain Monte Carlo (MCMC). Modern MCMC algorithms enable exploration of high-dimensional posterior distributions, have robust failure diagnostics (Vehtari et al. 2021) and can incorporate fast numerical solvers, thereby making inference feasible for Bayesian kinetic models. Nonetheless, the kinetic modelling literature reports an aversion to MCMC, rooted mainly in concerns about sampling time and the presumed difficulty of implementing the required statistical model (Andreas Raue et al. 2013; Saa and Nielsen 2016). We are only aware of two recent attempts to implement a Bayesian kinetic modelling approach using MCMC. Stapor et al. (2018) fitted detailed kinetic models using relatively inefficient MCMC algorithms that do not scale well to high dimensional parameter spaces limiting the scope of modelling. Conversely, St. John et al. (2019) utilises an efficient sampling algorithm but uses approximate kinetics, namely lin-log kinetics (Visser and Heijnen 2003), limiting the scope of interpreting parameters and inferring cellular behaviour in experimental conditions outside the reference dataset.
There have also been efforts to implement Bayesian inference for kinetic models without the use of MCMC. Examples of alternative inference methods include variational inference (St. John et al. 2019), rejection sampling and approximate Bayesian computation (Saa and Nielsen 2016) and Laplace approximation, in which the Fisher information matrix is used to calculate a normal approximation around the maximum a posteriori parameter configuration (Liebermeister and Noor 2021; Gopalakrishnan, Dash, and Maranas 2020; Stapor et al. 2018; Andreas Raue et al. 2013). Non-MCMC-based Bayesian kinetic models suffer from a lack of reliable diagnostic tools for verifying that their results approximate the target posterior distribution. This is a problem because realistic kinetic models tend to induce highly correlated, non-Gaussian, joint probability distributions (Gutenkunst et al. 2007; Stapor et al. 2018).
Our application Maud provides a Bayesian kinetic model that combines biologically realistic mechanistic accuracy—including accurate rate laws, post- translational modification and thermodynamics—with fast, robust MCMC sampling using adaptive Hamiltonian Monte Carlo. Further, Maud is a general-purpose application that can be used to fit a wide range of Bayesian kinetic models.
2 Results and Discussion
To demonstrate our application’s capabilities, we used Maud to analyse an artificial dataset based on the human methionine cycle. We generated this dataset using Maud, by simulating measurements for a set of training and validation conditions based on plausible parameter values. Next, we performed posterior inference for the training measurements and prediction for both training and validation measurements.
We investigated Maud’s sensitivity to missing measurements by comparing the results of fitting a full dataset with an intentionally incomplete dataset. To demonstrate why a full Bayesian approach is preferable to an approach based on a Laplace approximation of the posterior distribution, we compared the results of analysing a representative metabolic network using both methods.
Finally, we dug into our results to find out what our full dataset methionine model learned about the contributions of different regulatory factors to the flux through GNMT, an important reaction. This analysis illustrates how Maud can be used to generate actionable insights about metabolism without the need for further statistical analysis.
The methionine cycle, illustrated in Figure 1, is a fundamental pathway in human metabolism, whose intermediate metabolites participate in a variety of mechanisms which must compete for the same resources. Due to this competition, as well as the fact that all the functions occur simultaneously, the methionine cycle is highly regulated, with 6 known allosteric effectors (Martinov et al. 2000; Nijhout et al. 2006; Korendyaseva et al. 2008). This complex regulation means that quantitative modelling of the methionine cycle requires a detailed kinetic model: this is why we chose it as a case study for Maud.
2.1 Dataset and model specification
The simulated dataset and underlying kinetic model that we used for our analysis can be found at https://github.com/biosustain/Methionine_model/tree/main/data/methionine and is described in supporting information section 3.
We constructed a kinetic model of the methionine cycle in Maud’s format using the description in Korendyaseva et al. (2008). The ordinary differential equation system describing this model is
\[\begin{align} \frac{d[met]}{dt} &= v_{Influx} - v_{PROT} - v_{MAT} +v_{MS} + v_{BHMT} \label{eq-meth-ode} \\ \frac{d[amet]}{dt} &= v_{MAT} - v_{GNMT} - v_{METH} \nonumber \\ \frac{d[ahyc]}{dt} &= v_{GNMT} + v_{METH} - v_{AHC} \nonumber \\ \frac{d[hyc-L]}{dt} &= v_{AHC} - v_{CBS} - v_{MS} - v_{BHMT} \nonumber \\ \frac{d[5mthf]}{dt} &= v_{MTHFR} - v_{MS} \nonumber \end{align}\]
After specifying the qualitative aspects of the kinetic model, we selected parameter values to use as ground truth by Monte Carlo sampling using a previous model of the methionine cycle as a starting point (see Saa and Nielsen (2016) for this model).
We used these parameters to simulate steady states in a range of plausible experimental conditions, again using Saa and Nielsen (2016) as a starting point. These steady states were then used to generate simulated measurements using the measurement model.
For enzyme and metabolite concentration measurements we specified a standard deviation of 0.1 on natural logarithmic scale, corresponding to approximately 10% measurement error. For each reaction measurement a measurement standard deviation of approximately 10% of the simulated value.
These measurement error specifications are somewhat optimistic considering the many sources of variation and uncertainty affecting quantitative proteomics, metabolomics and fluxomics analyses, but are a reasonable first approximation to a realistic set of measurements.
For our main model run, we assumed that all metabolite and enzyme concentrations were measured, and that there was a reaction measurement for each of the network’s free fluxes.
The simulated experiments and measurements were split into a training and a validation dataset in a way that achieved a large difference in flux between the two categories. This was done to evaluate whether the fitted model is able to extrapolate to conditions well outside the training dataset rather than merely interpolating between the training data without necessarily learning the system.
We constructed inputs in Maud’s format for each of the analysed datasets, based on the scenario that the true kinetic model was known except for parameter values, which needed to be inferred from the training data and priors. These inputs can be found at https://github.com/biosustain/Methionine_model/tree/main/data.
2.2 Posterior inference
The prior distributions and corresponding true parameter values used in our case study are shown in supporting information section 3.2. They were chosen to reflect a plausible pre-experimental information state. In 7 cases, the marginal prior distribution for a parameter disagrees with the true parameter values used to generate the data; a similar situation is likely to occur in practice due to in vivo vs in vitro measurement differences.
Running standard diagnostic checks indicated that the samples we generated were from the target posterior distribution. The improved \(\hat{R}\) statistic (Vehtari et al. 2021) for every variable of interest was within 2% of 1, indicating appropriate mixing within and between Markov chains. Additionally, the number of effective samples was high, indicating that we generated enough posterior samples to support inferences about the bulks of the distributions of the sampled parameters. Furthermore, we observed no post warm-up divergent transitions, indicating that the sampler was able to transform the log-posterior distribution, avoiding any regions with excessive curvature that might inhibit exploration via HMC.
Posterior predictive checking indicated that our model achieved a good fit to the simulated reaction and metabolite concentration measurements, as shown by the graphs in the top row of Figure 2.
Analysis of the posterior distributions for kinetic parameters indicated that these are highly correlated. The marginal posterior distributions for most kinetic parameters did not shrink significantly compared to the corresponding marginal prior distributions, even though these parameters’ joint posterior distribution contained enough information to make accurate out of sample predictions. In some cases, there were two-dimensional correlations such as the one shown in the bottom left of Figure 2; in this case the marginal distribution of the two parameters is roughly banana shaped. More commonly, however, two-dimensional pair plots were insufficient to reveal the underlying correlation structure, as seen in the bottom-right plot in figure Figure 2, which shows two marginally independent parameters.
These results show that Maud can achieve an adequate fit to a realistic pathway-sized dataset. This was achieved without fixing the marginal values of kinetic parameters: the information required to make good predictions was contained in the correlation structure of the joint posterior distribution. This finding is consistent with previous analyses of biological systems that found they are “sloppy”, that is, sensitive to parameter combinations rather than marginal parameter values, with important combinations, scales and regions of sensitivity being difficult to ascertain in advance (Gutenkunst et al. 2007; Poirier 1998).
The question naturally arises whether the crucial high-dimensional parameter correlations are linear or non-linear. This is relevant to the question of model performance, as linear correlations are easier to correct for. A linearly correlated posterior space would also be easier to summarise. We address this question in the next section.
2.2.1 Comparison with Laplace approximation
This current case study illustrates the type of kinetic model and dataset that Maud can fit. The model we analysed has 10 reactions, 5 state variables and 212 parameters. Generalising from our ability to fit this model in a reasonable time using Maud, we expect that Maud can be used to fit realistic Bayesian models of approximately the same size, but not, for example, genome-scale kinetic models. To fit larger models, faster steady state solving methods or alternative inference algorithms will be required. The Laplace approximation, in which the Fisher information matrix is used to calculate a normal approximation around the maximum a posteriori parameter configuration, is a popular strategy for generating approximate posterior samples while avoiding a full MCMC approach.
Using the recently implemented Laplace approximation in Stan, we were unable to generate approximate posterior samples for our main methionine cycle case study using Laplace approximation, as the algorithm could not recover from solver failures caused by unrealistic parameter configurations. This is not an inherent issue of the Laplace transformation but highlights some of the challenges around solving \(\eqref{eq-3}\) for real problems.
In order still to assess the potential use of the Laplace approximation, we made a comparison for a simpler model. This input can be found at https://github.com/biosustain/Methionine_model/tree/main/data/example_ode. MCMC sampling for this simpler model yielded 800 samples in 625 minutes, while Laplace sampling yielded the same number of samples in only one minute. The diagnostics indicated that our algorithm was able to find the maximum a posteriori parameter configuration, approximate the Hessian and use these quantities to generate approximate posterior samples. The results can be found at https://github.com/biosustain/Methionine_model/tree/main/results.
Comparing the samples generated using each method shows that the Laplace approximation does not provide a good approximation to the true posterior distribution in this case (Figure 3). As can be seen from the top left plot, the total log probability density is clearly different and this was confirmed by a test of the equality of two empirical univariate distributions (Kolmogorov-Smirnov, \(p=1.7\times10^{-65}\)).
The difference between the Laplace approximation output and the true posterior distribution manifests not only in the parameter space, but also in the measurement space for both training data (upper right) and validation data (lower left). The lower log likelihood values indication that the modelled and measured values are further away and that the Laplace approximation yielded significantly worse predictions than the true posterior, even for the training data.
To further explore why this is the case we compared samples from the true posterior and the Laplace approximation for the pairwise marginal distributions of two Michaelis-Menten constants \(K_{m}^{A,r1}\) and \(K_{cat}^{r1}\) (Figure 3, bottom right). This comparison demonstrates that the Laplace method is not able to capture the correct relationships between parameters’ distributions.
This result shows that MCMC, while slower than Laplace approximation, is required for posterior inference in this case. We expect that the Laplace method will produce worse approximations the more complex the target model. Since the approximation is already unacceptable for our simple test model, we recommend that Maud users use MCMC sampling in preference to Laplace approximation if possible when fitting realistic Bayesian kinetic models.
Our results here also provide circumstantial evidence that the parameter correlations in Bayesian kinetic model posteriors tend to be non-linear, as a posterior with only linear correlations would likely be more germane to Laplace approximation. A conclusion that we drew from this analysis was that the results of fitting our model cannot be summarised simply, for example by fitting a multivariate normal distribution to the posterior draws. We therefore recommend that Maud users store the full set of MCMC draws rather than using such an approximation. This does not preclude the possibility that there is an alternative, more compact, way to summarise the results of Bayesian kinetic model inference; we leave research into this topic to future work.
2.3 Effect of missing metabolite concentration measurements
To gain insight into our model’s robustness to missing measurements, we also performed a model run with the same 6 experimental datasets, but with measurements of the metabolite S-Adenosyl-L-homocysteine, or “ahcys” removed. Since ahcys regulates three enzymes in the methionine cycle, including one enzyme which is also thermodynamically regulated, we expected the removal of these measurements to yield interesting results.
Comparing model runs with and without the ahcys measurements showed that Maud can produce sensible results even from incomplete metabolomics data.
As might be expected, the model with missing measurements did not correctly infer the missing ahcys concentrations (Figure 4). Nonetheless, the remaining measured metabolites were still well predicted, suggesting that information about the network is still preserved despite the missing measurements. Comparison of flux measurements in both models also indicated that removing the ahcys measurement did not result in catastrophic model failure.
The missing measurements did affect Maud’s ability to infer parameter values correctly (Figure 4, lower left). The model with full metabolomics learned the true value for the displayed dissociation constant, despite this value being far from the mean of the corresponding marginal prior distribution. In contrast, the model with missing measurements stayed in the neighbourhood of the prior.
This result is reassuring because not having access to all measurements is a common situation in multi-omics studies. For instance, measuring all metabolites in a pathway can be infeasible because of limitations of mass spectrometers, availability of standards, column effects, and compartmentalisation. However, provided that sufficient information is available from other sources, our approach can produce sensible results from incomplete metabolomics data.
2.4 Application to regulatory understanding
To demonstrate how Maud’s output can be used to yield useful metabolic insights we used the results of our case study to explain why the flux of the enzyme GNMT is higher in dataset 1 than in dataset 12 (Figure 5). GNMT is an irreversible enzyme that is homotropically activated by its substrate, competitively inhibited by its product and heterotropically inhibited by 5,10- methylenetetrahydrofolate (mlthf). The complex regulation makes it the ideal test case to elucidate regulatory changes.
Regulation can be separated into enzyme abundance, allostery and saturation, and we can plot the marginal posterior distribution of the ratio of the corresponding regulatory component in dataset 1 compared with dataset 12 (Figure 5, right panel). A positive value indicates that the component was increased in dataset 1 relative to dataset 12, with 0 indicating no difference. The probability, according to our model, of the component acting in each direction is given by the relative area under the curve on each side of the zero point.
Our model correctly inferred that saturation and allosteric effects were the main drivers of regulation between the two datasets in this case, with the curves for each component aligning with the ground truth shown in red.
Importantly, this form of analysis takes into account all modelled sources of uncertainty, including uncertainty about the measured values of the flux in each dataset. Our result shows that Maud could be used in this realistic case not only to explain an observed difference in fluxes, but also to provide a reasonable judgement as to the explanation’s robustness.
3 Methods
Maud is a command line application implementing Bayesian inference for a wide range of realistic kinetic models. Maud is written in Python (Van Rossum and Drake 2009), designed for use on Windows, macOS and Linux, is pip installable from the Python Package Index as maud-metabolic-models, documented at https://maud-metabolic-models.readthedocs.io and actively developed and maintained at https://github.com/biosustain/Maud/.
To use Maud, a user must first collate appropriate input information, represent it in files with Maud’s required formats (see supporting information section on input format). Maud’s command line interface provides commands for inference via a range of algorithms including adaptive Hamiltonian Monte Carlo, as well as commands for simulation and making out-of-sample predictions. Results are stored in files, using a structured, interoperable format.
As well as parameter values, Maud also performs inference for derived quantities including the components of the regulatory decomposition described below in \(\eqref{eq-decomposition}\), log likelihoods, simulated measurement values and metabolic control analysis coefficients (Kacser and Burns 1973).
In the rest of this section we describe Maud’s kinetic model at high level and discuss Maud’s statistical model, implementation and how it solves the crucial steady state problem. Further details about Maud’s kinetic model can be found in the supporting information.
3.1 Kinetic model
Maud’s kinetic model decomposes into factors:
\[\begin{equation} F(C;\theta) = Enzyme\cdot k_{cat}\cdot Reversibility \cdot Saturation \cdot Allostery \label{eq-decomposition} \end{equation}\]
Each of the terms on the right-hand side of \(\eqref{eq-decomposition}\) is a function of \(C\) and \(\theta\). This idea is taken from Noor et al. (2013). The terms usefully gather physically meaningful and conceptually distinct factors contributing to reaction fluxes. \(Enzyme\) captures the effect of enzyme concentration, \(k_{cat}\) that of enzyme efficiency, \(Reversibility\) quantifies thermodynamic effects, \(Saturation\) the effect of enzyme availability and \(Allostery\) the effect of post translational modifications.
We used the model of enzyme saturation from Liebermeister, Uhlendorf, and Klipp (2010) and the generalised Monod-Wyman-Changeux model of Allosteric regulation introduced in (Monod, Wyman, and Changeux 1965; Changeux 2013; Popova and Sel’kov 1975, 1979) and used more recently in Matos et al. (2022). To capture the effect on enzyme activity of coupled phosphorylation and dephosphorylation processes we developed a new mathematical model inspired by the generalised MWC model of allosteric regulation. Full details of all mathematical aspects of Maud’s kinetic model can be found in supporting information section 2.
3.2 Statistical model
Maud uses Markov Chain Monte Carlo to perform posterior inference on a Bayesian statistical model. This section introduces these topics and then describes exactly what kind of Bayesian statistical model Maud uses.
3.2.1 Bayesian inference
Bayesian statistical inference analyses data by constructing a mathematical model with the following three components:
A measurement model or “likelihood” that probabilistically describes how likely any possible measurement set is given the true values of the measured quantities, i.e. a probability density function \(l: \mathcal{Y}\times\hat{\mathcal{Y}} \rightarrow [0, 1]\) where \(\mathcal{Y}\) is the set of all possible measurement sets and \(\hat{\mathcal{Y}}\) is the set of all possible true measurable values.
A deterministic process model that describes the measured quantities in terms of unknown, possibly unobserved parameters \(\theta\), i.e. a function \(d: \mathcal{\theta} -> \hat{\mathcal{Y}}\).
A prior model that probabilistically describes how likely any possible set of parameter values is, without considering any information included in the measurement model, i.e. a probability density function \(\pi : \mathcal{\theta} \rightarrow [0, 1]\).
Together these components determine a joint probability function that encapsulates the Bayesian statistical model, i.e. a function \(p: \mathcal{Y}\times\mathcal{\theta}\rightarrow[0,1]\) such that for any \(\theta\) and \(y\), \(p(\theta, y) = \pi(\theta)l(y\mid\theta)\). Substantive questions about the implications of the assumptions implicit in the model components can be formulated in terms of this joint density function.
In particular, questions about what the model implies, given a particular measurement set \(y\) can be formulated in terms of the conditional probability density function \(p_{\mid y}:\mathcal{theta}\rightarrow[0,1]\). Bayes’s theorem guarantees that this conditional density is proportional to the product of the prior and likelihood, i.e. \(p_{\mid y}(\theta) \propto \pi(\theta)l(y\mid\theta)\) See (Gelman, Carlin, et al. 2020, Ch 1) for an extended discussion of the theory of Bayesian statistical inference.
The principal computational challenge for Bayesian statistical inference is to approximate the values of integrals of the joint density function conditional on a known measurement assignment \(y\), also known as the posterior distribution. For most non-trivial statistical models the posterior distribution cannot be integrated analytically, so numerical approximation is required. Markov Chain Monte Carlo or “MCMC” is a popular method for addressing this problem that uses a Markov chain with known properties to generate samples from a target probability distribution that can be used to perform Monte Carlo integration. Hamiltonian Monte Carlo makes it possible, given an appropriate choice of hyperparameters, to efficiently generate samples from even a high dimensional continuous posterior distribution, by calculating the gradients of the log-scale joint density function. Adaptive Hamiltonian Monte Carlo, as used by Stan, uses a range of well-tested methods to tune these hyperparameters, allowing efficient MCMC for high dimensional posterior distributions with minimal user input.
Alternatives to MCMC for numerically approximating integrals over posterior distributions include grid sampling (Gelman, Carlin, et al. 2020, Ch. 10), rejection sampling (Gelman, Carlin, et al. 2020, Ch. 10), importance sampling (Vehtari et al. 2022) and distributional approaches including variational inference and Laplace approximation (Gelman, Carlin, et al. 2020, Ch. 13).
The Laplace approximation is of particular interest because, as mentioned above, this has been used for approximate Bayesian inference in a similar context to Maud. Laplace approximation of a posterior distribution works by first finding the mode of the posterior distribution, i.e. the ‘maximum a posteriori’ parameter configuration with the highest posterior density. Next, the second derivative of the posterior distribution at this point is found and used to generate a normal distribution that either approximates the target distribution, or can in turn be used to generate such an approximation.
Maud employs adaptive Hamiltonian Monte Carlo to perform posterior inference for a Bayesian statistical model where the process model is the kinetic model described above. The next two subsections describe Maud’s prior and measurement models.
3.2.2 Prior model
Maud’s prior model includes unknown parameters corresponding to quantities in the kinetic model that are assumed to be unknown, other than steady state metabolite concentrations and fluxes, which are derived from the values of other parameters by solving the steady state problem. See table 1 in this paper’s supporting information for a description of all these parameters and their dimensions.
Except for metabolites’ standard condition Gibbs energy changes of formation, Maud uses independent normal prior distributions for parameters that can in principle be both negative and positive. For parameters that are constrained to be positive, Maud uses independent log-normal distributions. Formation energy parameters have a multivariate normal prior distribution. Location, scale and covariance parameters for all these prior distributions can be selected freely by the user.
3.2.3 Measurement model
Maud’s measurement model considers three types of measurement: metabolite concentration measurements, enzyme concentration measurements and flux measurements, represented by vectors \(𝑦^{𝑐𝑜𝑛𝑐}\) \(𝑦^{𝑒𝑛𝑧}\) and \(𝑦^{𝑓𝑙𝑢𝑥}\) respectively.
All measurements are specific to an experimental condition; that is, a case where the true state of the network, including knockouts, boundary conditions and state variables as well as kinetic and thermodynamic parameters, can safely be assumed to be the same. Maud’s statistical model allows for arbitrarily many experimental conditions, and for any measurable quantity to be measured any number of times in any condition.
Metabolite and enzyme measurements are intended to represent the results of quantitative metabolomics and proteomics experiments. The likelihood functions for such measurements are:
\[\begin{align} y_i^{conc} &\sim LN(\ln{\hat{y}_i^{conc}}, \sigma_i^{conc})\label{eq-yconc} \\ y_i^{enz} &\sim LN(\ln{\hat{y}_i^{enz}}, \sigma_i^{enz})\label{eq-yenz} \end{align}\]
Both equations are log-normal generalised linear models with a standard link function (the natural logarithm \(\ln\)) and known standard deviation $ _{𝑐𝑜𝑛𝑐}$. The use of this measurement model is motivated by the consideration that concentrations are constrained to be non-negative, so the measurement model should avoid assigning positive probability mass to negative metabolite concentration values. In addition, we expect the precision of most metabolomics and proteomics experiments to be roughly proportional to the true value of the measured quantity, which supports a measurement model with constant coefficient of variation. The measurement standard deviations \(\sigma_{𝑐𝑜𝑛𝑐}\) and \(\sigma_{𝑒𝑛𝑧}\) are assumed to be known exactly for simplicity; plausible values can be elicited by considering the likely coefficient of variation of the measuring apparatus. Our measurement model improves on analyses of metabolomics and proteomics data that assume a regression model with normally distributed errors, whether explicitly using a standard linear model or implicitly using ordinary least squares fitting.
Flux measurements, representing the results of quantitative fluxomics analyses, are modelled using a likelihood function from a standard linear regression model \(\eqref{eq-yflux}\). Flux measurements can be obtained from isotope labelling experiments using metabolic flux analysis, for example as described in (Young 2014). When entering flux measurements, it is important only to specify measurements for a network’s free fluxes, as the values of some steady state fluxes in a metabolic network are constrained by others, with the result that dependent fluxes cannot typically be measured separately. If measurements of multiple dependent fluxes are entered, information will inappropriately be double counted.
\[\begin{equation} y_i^{flux} \sim LN(\ln{\hat{y}_i^{flux}}, \sigma_i^{flux})\label{eq-yflux} \end{equation}\]
The use of independent measurement models for metabolite, enzyme and flux measurements carries an implicit assumption that there are no systematic correlations in the measurement errors. This choice was motivated by simplicity - it would be better to use a model with potentially correlated measurements. Similarly, it would be preferable to include measurement errors as model parameters, thereby avoiding possible bias due to incorrect assessments of measurement accuracy. However, we chose to use a simpler measurement model to avoid the complexity and potential fitting issues that these changes would entail.
Finally, the reader may wonder why Maud uses a linear regression model for reaction flux measurements even though this creates the potential for erroneous double counting and requires non-trivial upstream modelling, as intracellular fluxes typically cannot be measured directly. Instead, fluxes are typically inferred from isotope labelling experiments using metabolic flux analysis: see Dai and Locasale (2017) for more about this method. Ideally Maud’s measurement model for fluxes would extend from fluxes to the results of potential labelling experiments, thereby removing the need for upstream analysis and avoiding any double counting. This option has not yet been pursued, again for the sake of simplicity.
3.3 Implementation
Maud uses the Python library click (Click Developers 2022) to implement a command line interface. The command line interface loads input files as Python dictionaries, which are parsed using the Python library toml (Pearson 2020) and then validated and converted into structured MaudInput objects using Pydantic (Pydantic developers 2022). Maud’s statistical model is implemented in the probabilistic programming language Stan (Carpenter et al. 2017) and accessed using the interface cmdstanpy (Stan Development Team 2022). For posterior sampling, Maud uses the MaudInput to create an input file for Stan and then trigger posterior sampling using adaptive Hamiltonian Monte Carlo. See Betancourt (2018) for more about this algorithm.
When sampling is complete, Maud converts to the output into the standard format InferenceData using the Python library arviz (Kumar et al. 2019) and saves it as a json file, along with some information for debugging. This format allows for easy checking of MCMC diagnostics including divergent transitions, effective sample size and the improved \(\hat{R}\) statistic proposed in Vehtari et al. (2021).
3.4 Solving the steady state problem
Maud’s parameters are connected with measurable quantities via the steady state equation \(\eqref{eq-3}\). Posterior sampling using adaptive Hamiltonian Monte Carlo requires repeatedly solving \(\eqref{eq-3}\) and finding the gradients of its solution with respect to all model parameters. This must be done numerically as analytic solutions are not available for realistic kinetic models.
The speed at which this problem can be solved is tightly coupled with the size and complexity of metabolic network that can practically be modelled. See Timonen et al. (2022) for more about considerations involved in this kind of modelling.
To solve \(\eqref{eq-3}\) and find its gradients, Maud uses a hybrid method involving two numerical solvers from the SUNDIALS suite (Serban and Hindmarsh 2005): CVODES and IDAS via their interface from Stan. The hybrid method follows that proposed by Margossian (2018), and involves numerically evolving the ODE system for a short period of time, then using the difference between the evolved and starting concentrations as the target for a numerical algebra solver.
4 Supporting information
The analyses described in this paper, as well as instructions for reproducing our results, figures and paper, can be found at https://biosustain/Methionine_model.
The document supporting_information.pdf contains the following information:
- Description of Maud’s input format.
- Description of Maud’s kinetic model including rate equations as well as parameters and their dimensions.
- Details of the procedure used to generate our case study results, including generating the artificial case study datasets, specifying priors and carrying out computation.
- Specification of the prior distribution used for the case studies.
5 Acknowledgements
This research was funded by the Novo Nordisk Foundation (Grant numbers NNF20CC0035580 NNF14OC0009473).